#!/usr/bin/perl -w
use strict;
use List::Util qw[min max sum];
use Class::Struct;

## Input1: precursor m/z and charge states
## Input2: tandem MS m/z and peak intensity
## Basic idea: the fragments are generated from all possible structres fit the precursor m/z and charge states.

## INPUT in a batch mode.
## The basic format is:
## PRE_M/Z CHARGE
## PRO_MS  ABD

my $dirname = "";
opendir(DIR, $dirname) or die "Cannot open the directive $dirname:$!";
my @files = grep(/\.dta$/,readdir(DIR));
closedir DIR;

## Element information ##
## Monoisotopic mass   ##

###################################
## Parse the deisotoped dta file ##  
###################################

## Get composition for the precursor ion ##
#my @compo_group = 0 x 5; # Composition: (delta-HexA, HexA, GlcN, SO3, Ac)



## Generate all possible fragments based on compositions ##
## 1. All possible structures ##
## 2. All fragments deduced from those structures ##
## 3. Match the peak list against the theoretical fragments ##
## Tolerance should be taken into consideration ##

################################
## 1. All possible structures ##
################################
## Assume the composition is (1,2,3,2,1) and the charge state is 2
our @compo = (1,2,3,2,2);
## dp number is defined by the first 3 digits.
our $dp = &sum(@compo[0..2]);
## In order to represent the complicate structure, a new nomenclature can be used.
## HexA: 2C(S)
## GlcN: 2C(Ac/S),3C(Ac),6C----Ac is defnied as digit 3, all others are defined as the position number. 
my @modi_chain;
## HexA: 0 | 2
## GlcN: 0 | 1 | 2 | 3 | 6 

## Define the position of acetation: select $compo_group[4] from $dp/2
#our @hexa_val = (0,1);
#our @glcn_val = (1,10,100,2);
our @glycan_val = (
		   ['0','2'],
		   ['0','1','2','3','6','13','16','23','26','36','136','236'] # 1 -- Ac modified.
		);
my $candidate_chain; # Reference to all candidate chains
my $pos = 0;     

## Record the total number of sulfate/acetate
my @sum_info = (0,0);

## Update candidate chain with possible positions.
$candidate_chain = ExpModi(0,\@modi_chain,$candidate_chain,\@sum_info);
open(OT,">possible_structure.txt") or die "Cannot write to this file:$!";
foreach my $chain (@{$candidate_chain}){
	my $string = join "\t", @{$chain};
	print OT "$string\n";
}
print "Done!\n";

close(OT);

####################################################
## 2. All fragments deduced from those structures ##
####################################################

## 2.1 Calculate the mass of several chemical groups
my %MonoTable = (
		'C' => ,
		'H' => ,
		'O' => ,
		'S' => ,
		'N' => 
		);
# HexA: C6H10O7
my $M_hexa = 6 * $MonoTable{'C'} + 10 * $MonoTable{'H'} + 7 * $MonoTable{'O'};
# GlcN: C6H13NO5 #
my $M_glcn = 6 * $MonoTable{'C'} + 13 * $MonoTable{'H'} + 5 * $MonoTable{'O'} + $MonoTable{'N'};
# Water: H2O #
my $M_h2o = 2 * $MonoTable{'H'} + $MonoTable{'O'};
# Sulfate group: SO3 #
my $M_sg = $MonoTable{'S'} + 3 * $MonoTable{'O'};
# Acetate group: C2H2O #
my $M_ac = 2 * $MonoTable{'C'} + 2 * $MonoTable{'H'} + $MonoTable{'O'};
## Connecting Oxygen.
my $M_o = $MonoTable{'O'};

## Calculate the mass of structural units
# HexA
my ($N1_c2h2o2, $N1_ch2o, $N1_o) = (1, 4, 1);
# GlcN
my ($N2_c2h4o, $N2_ch2o, $N2_ch3n, $N2_o) = (1, 3, 1, 1);

my @links = (1,4);
my @all_types;

## Only intested in the overall mass. The information might be known from the precursor information. ##
## All the modification mass will be deducted from the overall mass.
my $candidate_mass = int($dp/2) * $M_hex + ($dp-int($dp/2)) * $M_glcn;

## Backbone mass.
my $mass_loss=0;
my @target_chain = (1,(0) x (2 * $dp -1),-1);
my @chain_mass = (0) x 2*$dp+1;
my @chain_component;

for(my $pos=0; $pos<@target_chain; $pos++){
	## Do nothing with the terminal: mass is zero.
	next if $pos==0 || $pos==$#target_chain;
	## The connecting oxygen atom.
	if($pos % 2 == 0){
		$chain_mass[$pos] = $M_o;
		next;
	}
	if($pos == 1){    # Unsaturated HexA.
		$chain_mass[$pos] = $M_hexa - $M_h2o - $M_oh;
	}
	if(($pos-1)%4 == 0){  # HexA
		$chain_mass[$pos] = $M_hexa - $M_h2o;
	} else {	      # GlcN
		$chain_mass[$pos] = $M_glcn - $M_h2o
	}
		
}

## Modification of chain and corresponding mass change.
our @chain_lib;
my $index=0;
foreach my $chain (@{$candidate_chain}) {
	## Create a temporary variable to operate the mass.
	my @new_chain = @chain_mass;
	##\\ The exact value will be filled later.
	my @hex_ring = ();
	my @glc_ring = ();
	## For each chain, a variety of modification will be applied, combined with sulfate/water loss.
	for(my $i=0; $i<@{$chain}; $i++){
		##\\ Convert the string to number
		my @modi = split '', $chain->[$i];
		## If there is no modification, go to the next number.
		if($modi[0]==0){    # Please take care of the digit or string.
			next;	
		}
		## If there is some modification...
		for(my $j=0; $j<@modi; $j++){    # For each digit cluster.
			
			if($modi[$j] == 1){ # Acetation on the 2-N. There is no need to consider the ring number.
				## Add Ac mass and reduce water mass.
				## Update the mass.
					
				$chain_mass[$i] += $M_ac;
				$chain_component[$i]->[$j] += $M_ac; 
			} else {
				## Add Sg mass and reduce water mass.
				$chain_mass[$j] += $M_sg;	
				$chain_component[$i]->[$j] += $M_sg;
			}
		}	
		
	}
	$chain_lib[$index++] = cleavage_edge($chain_component[$i]);
}

## Test the results of step 2.## 
open(FL,">fragments.txt") or die "Cannot open file:$!";
my $chain_order = 0;
foreach my $chain (@{candidate_chain}){
	
	my $string = join "\t", @{$chain};
	print FL "::$string\n";	
	
	foreach(@{$chain_lib[$chain_order++]}){
		## Think twice of the exact format of deep_sum.	
		my $sum_ms = deep_sum(@$_);
		print FL "$sum_ms\t@$_\n";
			
	}
	
}
close(FL);

##############################################################
## 3. Match the peak list against the theoretical fragments ##
##############################################################
## 3.1 Read the peaklist ##
## m/z and z, deisotoped.##
my @peaks; 

## 3.2 Match m/z of peaks against the fragments.
## 3.2.1 Sort the masses of fragments for separate precursor ions.






## 3.2.2 Store the matched mass, calculate the hit number and percentage.


## 3.3 Evaluate the possibility of different fragments from precursor-wise comparison.
## 3.3.1 Compare the hit fragments, see the similarity across different precursors, especially the one with similar
## structure. Do some statisticla analysis using the data.
## I would expect the true structure show significantly high scores among the comparison.
## Score to discover and summarize the similarity.
## 3.3.2  

## Act as a filter to filt out impossible structures          ##k
## The control parameter should be passed into the subroutine ##
## to keep the intact of the function 			      ## 
sub BioSyn {

}

## A recursive subroutine to generate all possible modification combinations ##
sub ExpModi {
	## $index starts from 0 till $dp-1.
	my ($index,$modi,$candi,$info)= @_;
	my $type = $index % 2;  # 0 -- HexA and 1 -- GlcN
	
	## Go through all possible modification types and count the corresponding modification(s).
	foreach my $pattern (@{$glycan_val[$type]}){
		## Start of a new chain.
		$modi = [(0) x $dp] if $index==0;

		my $digit_count = length($pattern);
		my @init = split '', $pattern;
		my @temp = (0,0);
		
		if($init[0] eq '1'){ # Acetation.
			
			next if $info->[1] == $compo[4];
			$temp[0] = $info->[0] + ($digit_count - 1);
			$temp[1] = $info->[1] + 1;
			
			next if $temp[0] > $compo[3]; # Discard this possibility.
		#	@{$info} = @temp;

		} elsif($init[0] ne '1' && $init[0] ne '0') { # Sulfation
			$temp[0] = $info->[0] + $digit_count;
			$temp[1] = $info->[1];
			
			next if $temp[0] > $compo[3];
		#	@{$info} = @temp;
		}

		## Update $modi
		$modi->[$index] = $pattern;     
		
		## End of chain
		if($index == $dp-1) {
			my @temp_seq = @{$modi};
			push @{$candi}, \@temp_seq;
			next; 
		}

		$candi = ExpModi($index+1,$modi,$candi,\@temp);
	}
	
	return $candi;
}

## generate all possible cleavage
sub cleavage_edge {

	my $chain_mass = shift @_;
	
	## The cleavage is based on edges, so the numbers will be redefined.
	## Take account the terminal (begin/end), there will be 2*dp edges ranging from 0 -> 2*dp-1
	for(my $nre=0; $nre<2*$dp-1; $nre++){

		for(my $re=$nre+1; $re<2*$dp; $re++){

			my @temp_chain_mass = @$chain_mass;
			
			## All the mass before $nre and after $ne will be cleared.
			## Update the mass of @temp_chain_mass.
			@temp_chain_mass = ((0) x (0 .. $nre), $temp_chain_mass[$nre..$re-1], (0) x ($re+1..2*$dp))

			## The mass between $nre and $re will be calculated considering the modifications.
			## 1. Calculate the full mass.
			#my $overall_mass = 0;
			#foreach(@{$lib}[$nre..$re-1]){
				## Add the mass of the unit.
			#	$overall_mass += sum(@$_);
			#}	
			
			## Update the mass on the notch.
			@temp_chain_mass = @{cleavage_ring(\@temp_chain_mass, $nre, $re)};
			

		}

	}
	return $chain_cleavage;
}

sub cleavage_ring {
	
	my ($cleavage_chain, $nre, $re) = shift @_;
	## 2. Deduce the missing part.
	## Treat NRE and RE separatly.
	my ($nre_frg,$re_frg);
	for(my $i=0; $i<4; $i++){    # The first index of the cleavage.
		
		for(my $j=$i+2; $j<=5; $j++){  # The second index of the cleavage. 			
			
			my $nre_type = decide_type();
			my $re_type = decide_type();
			next if $nre_type eq 0; # This is an internal cleavage.
			next if $re_type eq 0;	
			
			
			## Decide which part will be removed.
			## Treat NRE and RE differently.
			## a. NRE ##
			next if $nre==0 && $re == 2*$dp;

			## Update the mass based on cleavages.	
			if($links[1] > $i && $links[1] <= $j){
				@{$cleavage_chain->[$nre]}[$i..($j-1)] = (0) x @{$cleavage_chain->[$nre]}[$i..($j-1)];
				@{$cleavage_chain->[$re]}[0..$i] = (0) x @{$cleavage_chain->[$re]}[0..$i];
				@{$cleavage_chain->[$re]}[($j+1)..5] = (0) x @{$cleavage_chain->[$re]}[($j+1)..5];
				
			} elsif($links[0] > $i && $links[1] > $j) {
				@{$cleavage_chain->[$nre]}[0..$i] = (0) x @{$cleavage_chain->[$nre]}[0..$i];
				@{$cleavage_chain->[$nre]}[($j+1)..5] = (0) x @{$cleavage_chain->[$nre]}[($j+1)..5];
				@{$cleavage_chain->[$re]}[$i..($j-1)] = (0) x @{$cleavage_chain->[$re]}[$i..($j-1)];
			}
			push @{$nre_frg}, \@{$cleavage_chain->[$nre]};
			push @{$re_frg}, \@{$cleavage_chain->[$nre]};

		}
	}		
		
	return $cleavage_chain;	
	
}

sub decide_type {    # Params: $i, $j.
	## For terminal unit, any cleavage is X-ring.
	## For others, the X-ring depends on the cleavage sites and direction.
	my @index = @_;
	my $flag = 1;
	if($index[2] == 0 || $index[2] == 2*$dp-1) {
		return $flag;
	}
	
	if(($links[0]>$index[0] && $links[1]<=$j) || ($links[0] <=$i && $links[1] > $j)) {
		$flag = 0;
	}	
	return $flag;	
	
}
sub sum {
	my @list = @_;
	my $sum = 0;
	$sum += $_ foreach @list;
	return $sum;
}

